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q 1 Abstract. We present a new numerical scheme for solving the advection equation 

q | and its application to the Vlasov simulation. The scheme treats not only point val- 

^ . ues of a profile but also its zeroth to second order piecewise moments as dependent 

q 1 variables, and advances them on the basis of their governing equations. We have de- 

• ^ \ veloped one- and two-dimensional schemes and show that they provide quite accurate 

solutions compared to other existing schemes with the same memory usage. The two- 
dimensional scheme can solve the solid body rotation problem of a gaussian profile 
with little numerical diffusion. This is a very important property for Vlasov simulations 
'— — of magnetized plasma. The application of the scheme to the electromagnetic Vlasov 

simulation of collisionless shock waves is presented as a benchmark test. 

o. 

<N ■ 

^nI" \ 1. Introduction 

(N; 

The kinematics of collisionless plasma has been studied in a wide variety of fields, such 
as in laboratory plasma physics, space physics, and astrophysics. Evolution of collision- 
less plasma and self-consistent electromagnetic fields is fully described by the Vlasov- 
Maxwell (or Vlasov-Poisson) equations. Thanks to recent development in computa- 
^ . tional technology, self-consistent numerical simulations of collisionless plasma have 

been successfully performed from the Vlasov-Maxwell system of equations. 

One of numerical simulation methods for collisionless plasma is the so-called 
Vlasov simulation, in which the Vlasov equation is directly discretized on grid points in 



phase space. Compared to the most popular Particle-In-Cell (PIC) method (IBirdsall & Langdon 



19911) . the Vlasov simulation is free from the statistical noise inherent to the PIC 
method. This advantage can allow us to study in detail such as wave-particle inter- 
action, particle acceleration, and thermal transport processes, in which a high energy 
tail in the velocity distribution function plays an important role. On the other hand, 
the Vlasov simulation requires a highly accurate scheme for the advection equation in 
multidimensions, to preserve characteristics of the Vlasov equation (i.e., the Liouville 
theorem) as much as possible. It also requires larger computational cost than the PIC 
method. 
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A number of advection s chemes have been proposed for the application to the 
Vlasov simulation thus far (e .g. JCheng & KnorrHl976l : rNakamura & Yabdl999l : lFiibet et al 



2001; Mangen ey et al.ll2002h . Although the schemes have been succeeded especially in 
applying to the electrostatic Vlasov-Poisson simulation, the application to the electro- 
magnetic Vlasov simulation of magnetized plasma is still limited, mainly owing to the 
difficulty in solving the gyro motion around the magnetic field line. 

In this paper, we propose a new numerical scheme for the advection equation, 
specifically designed to solve the Vlasov equation in magnetized plasma. The scheme 
is briefly introduced in Section |2] Benchmark tests of the scheme and its application to 
the Vlasov simulation are presented in Section |3] Final ly, we summarize the pa per in 
Section 0] Details of the scheme have been presented in Minoshim a et all (1201 lh . 



2. Multi-Moment Advection scheme 

The present scheme considers the advection of a profile f(x, t) and its zeroth to second 
order moments denned as 



M m = ^- f x m fdx, (m = 0,1,2). (1) 
ml J 

In one dimension, their governing equations are written as 
df d 

17 + ^-("/) = 0, (2) 
at ox 

dM° C d 

— + J dx— (uf) = 0, (3) 

dM m i r d i r 1 

+ dx M yn/ = u^fdx, (m = 1,2), (4) 

at ml J ox (m - 1)! J 

where u is the velocity. Equations (0) and (O are exactly obtained by multiplying Equa- 
tion (O by x m /ml and then integrating over space. To solve a set of these equations, the 
one-dimensional scheme treats four dependent variables; the point value of the profile 
fi, and the piecewise moments, 

M il\/2 = -7 f *'"fdx, (m = 0, 1, 2) , (5) 



- f 

ml J Xi 



and constructs a piecewise interpolation for / in a cell with a fourth order polynomial, 
Fi(x) = ZjLi kCtj{x - Xj) k ~ l . The five coefficients Ctj are explicitly determined from 
the dependent variables at the upwind position as constraint. Then the variables are 
advanced on the basis of their governing equations ©-((U) with the semi-Lagrangian 
method. 

The two-dimensional scheme is designed in a similar way. It treats six dependent 
variables; the point value of the profile fij, and the piecewise moments in the x and y 
directions, 

M r + i/2j + i/2 = ^ I ' x m fdxdy, (m = 0,1,2), (6) 

hi 
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Figure 1 . One-dimensional linear advection of a gaussian profile. (a,b,c) Calcu- 
lation results with the MMA, the CIP-CSL2, and the PFC schemes (solid lines with 
symbols). Dashed lines are the exact solution. (d,e,f) Deviation of the calculation 
results from the exact solution. 



and constructs a piecewise interpolation for / in a cell with a quadratic polynomial, 
Fij(x,y) - Tq = \ ZLi lkCitjj(x - Xj) k ~ l (y - The nine coefficients Citjj- are ex- 

plicitly determined from the dependent variables at the upwind position as constraint. 
Th e scheme is termed as th e "Multi-Moment Advection (MMA)" scheme. For details, 
see Minoshima et all (|201lt) . 



3. Benchmark tests 



Figure Qla-c) shows the one-dimensio nal linear advect ion problem of a gaussian pro- 
file solved by the MMA, CIP-CSL2 (lYabe et all l200ll) . and PFC dFilbet et all [20011) 
schemes. Since the numbers of dependent variables are different among the three 
schemes (four for the MMA, two for the CIP-CSL2, and one for the PFC), we use 
different grid sizes so that the total memory usage is equal. The CFL number is 0.2. 
The MMA scheme (a) provides a quite accurate solution compared to other schemes 
(b,c). Figure QTd-f) shows the deviation of the calculation results from the exact solu- 
tion. The MMA scheme (d) is about fifty times better then other schemes (e,f). 

Figure I3a-f) shows the two-dimensional solid body rotation and advection prob- 
lem of a symmetric gaussian profile, 



d f < ^ d f , df 
¥ -0>-)>o)^ + (*-*o)^ 



0, / (x, y, t - 0) = exp 



(* - x f + (y - y ) 2 



2cr 2 



solve d by the MMA, CIP-CSL2 dTakizawa et al.l2002l) . and backsubstitution (|Schmitz & Grauer 
2006) schemes. This describes the rotation around (x,y) = (xo,yo), corresponding to 
the electric field drift motion for magnetized plasma. The parameters are (xo,yo,c) - 
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Figure 2. Two-dimensional solid body rotation and advection of a symmetric 
gaussian profile. (a,b,c) Calculation results after 50 rotations with the MMA, the 
CIP-CSL2, and the backsubstitution schemes. (d,e,f) Calculation results after 300 
rotations, (g) Temporal variation of the standard deviation cr. Solid, dashed, and dot- 
dashed lines are obtained from the MMA, CIP-CSL2, and backsubstitution schemes. 



(-0.05,0.1,0.1), and the simulation domain is [-0.5,0.5] in both directions. Since 
the numbers of dependent variables are different among the three schemes (six for the 
MMA, four for the CIP-CSL2, and one for the backsubstitution), we use the different 
numbers of grid points (34x34 for the MMA, 42x42 for the CIP-CSL2, and 84x84 for 
the backsubstitution) so that the total memory usage is equal. The time steps are 0.0047T 
for the MMA and CIP-CSL2 schemes, and 0.002tt for the backsubstitution scheme so 
that the CFL number is close among the three simulations. While other schemes show 
serious numerical diffusion after several tens of rotation periods, the MMA scheme 
completely preserves the profile after hundreds of rotation periods. Figure |2£g) shows 
the temporal variation of the standard deviation cr obtained by fitting the profile with 
the gaussian function. After 300 rotation periods, the standard deviation is increased 
by 0.1006 (MMA), 0.1272 (CIP-CSL2), and 0.1436 (backsubstitution). 

We apply the MMA scheme to the electromagnetic Vlasov-Maxwell simulation. 
The one-dimensional Vlasov-Maxwell system of equations is written as 

df s df s q s ( vxB\ df s 

-£+Vx^ 1 + —[E+ ~=0, (s = p,e), (7) 

ot ox m s \ c ov 



— = cVxB- Anj, — = -cVxE, /' = Jj q * \ vfsdv > 

' s=p,e 



(8) 



where E and B are the electric and magnetic fields, j is the current density, c is the speed 
of light, q s is the charge, m s is the mass, and the subscript s denotes particle species (p 
for protons and e for electrons). We assume the two dimensionality in velocity space, 
v = (v x , v y , 0), E = (E x , E y , 0), and B = (0, 0, B z ). The Vlasov equation is split into 
two equations in two-dimensional velocity and one-dimensional configuration spaces, 
which are advanced by the MMA and CIP -CSL2 schemes, r espectively. The Maxwell 
equation ([8]) is solved by the CIP scheme (IQgata et al.ll2006l). The tim e integration of 



the system is carried out in the same manner as lMinoshima et all (2011). 
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Figure 3. One-dimensional electromagnetic Vlasov simulation of perpendicular 
shock waves. (a,b) The electron phase space distributions in (x, v ex ) and (x, v ey ). 
(c,d) The proton phase space distributions in (x, v px ) and (x, v py ). (e,f,g) The elec- 
tromagnetic field (B z ,E x ,E y ) distributions. The velocity and the electric field are 
normalized by the bulk flow velocity and the motional electric field at the upstream, 
respectively. 



We simulate one-dimensio nal, strictly perpendicular collisionless shock waves 
(e.g.. iHoshino & Shimadall2002h . A high speed plasma is injected from the left-hand 
boundary and flows toward positive direction. The plasma carries the perpendicular 
magnetic field. At the right-hand boundary, the plasma is specularly reflected. As a 
result, a shock wave is formed and propagates in negative direction. Simulation pa- 
rameters are as follows; a proton to electron mass ratio m p /m e = 25, a ratio of the 
electron plasma to gyro frequency a) pe /a) ge = 100, electron and proton plasma beta 
values p e = p p = 1.0, and an Alfven Mach number of the upstream plasma flow is 
5.0. The velocity space domain is [-0.06c, 0.06c] for electrons and [-0.03c, 0.03c] for 
protons with 72 grid points in both the v x and v y directions. The configuration space 
domain is 20480/1/) with 1024 grid points (Ax = 20Ad), where Ad is the Debye length. 
The time step is At - 0.lna)~l. 

Figure [3] shows the electron phase space distribution ( J f e dv y , J f e dv x ), the proton 

phase space distribution ( J f p dv y ,J f p dv x ), and the electromagnetic fields (B z ,E x ,E y ) 
at oj ge t = 100k. An Alfven Mach number of the resulting shock wave is ~ 7.5 mea- 
sured in the shock rest frame. The simulation describes fundamental structures of the 
perpendicular collisionless shock. The plasma pressure and the magnetic field strength 
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rise at the shock front (x - 50), and subsequently oscillate due to the gyro motion of 
protons in the downstream region (x > 50). Around the front, the difference in inertia 
between electrons and protons produces the electrostatic potential (Figure Of), the so- 
called shock potential). Before the front, there is a gradual increase of the magnetic field 
strength (30 < x < 40), in which part of protons are reflected by the shock potential 
(Figure [He), the so-called reflected ions). We confirm that the Rankine-Hugoniot con- 
servation laws are satisfied at the shock. Even at this moment, the electron magnetic 
moment is well conserved in the downstream region, due to the fact that the MMA 
scheme can solve the solid body rotation with little numerical diffusion. 

4. Summary 

We have presented a new numerical scheme for solving the advection equation and the 
Vlasov equation. The present scheme solves not only point values of a profile but also 
its zeroth to second order piecewise moments as dependent variables, and advances 
them on the basis of their governing equations. We have developed one- and two- 
dimensional schemes, and have shown their high capabilities. The scheme provides 
quite accurate solutions compared to other existing schemes with the same memory 
usage. The two-dimensional scheme can solve the solid body rotation problem of a 
gaussian profile with little numerical diffusion. This is a very important property for 
Vlasov simulations of magnetized plasma. 

The application of the scheme to the electromagnetic Vlasov simulation of colli- 
sionless shock waves has been presented. In the simulation, we use co pe /u>g e = 100 and 
Ax = 20/1/). Although the grid size is much larger than the Debye length so that the 
Debye-scale structures can not be described, the simulation is stable and the meso-scale 
shock structures are well described. Since the grid size of an explicit PIC simulation 
is restricted to the Debye length, the PIC simulation requires large computational cost 
when ojpe/ojge is large. This is not the case for Vlasov simulations, unless Debye-scale 
structures are important. This advantage enables us to perform large-scale plasma ki- 
netic simulations with large u> pe /u> ge . 

Acknowledgments. We thank an anonymous referee for reviewing the manuscript. 
Part of Figures [T]and |2] are reproduced by permission of the Elsevier Inc. 
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